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Abstract: We discuss applications of time-delayed feedback control to delay-coupled neural 
systems and lasers, in the framework of the FitzHugh-Nagumo neuron model and the Lang- 
Kobayashi laser model, respectively. In the context of neural systems, we will point out some 
complex scenarios of synchronized in-phase or antiphase oscillations, bursting patterns, or 
amplitude death, induced by delayed coupling in combination with delayed self-feedback in 
simple network motifs. For optical systems, we will show that multiple time-delayed feedback, 
realized by a Fabry-Perot resonator coupled to the laser, provides a valuable tool for the 
suppression of unwanted intensity pulsations, and leads to stable continuous- wave operation. 
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Over the past decade control of unstable and chaotic 
states has evolved into a central issue in applied nonlinear 
science (Scholl and Schuster, 2008). Various methods of 
control have been developed since the ground-breaking 
work of Ott et al. (1990). One scheme where the control 
force is constructed from time-delayed signals (Pyragas, 
1992) has turned out to be very robust and universal to 
apply, and easy to implement experimentally. In this time- 
delayed feedback control the control signal is built from the 
difference s{t) — s{t — r) between the present and an earlier 
value of an appropriate system variable s. It is noninvasive 
since the control forces vanish if the target state (a periodic 
state of period r or a steady state) is reached. Thus the 
unstable states themselves of the uncontrolled system are 
not changed, but only their neighbourhood is adjusted 
such that neighbouring trajectories converge to it, i.e., the 
control forces act only if the system deviates from the state 
to be stabilized. 

This paper is organized as follows: In Sec. 1, we discuss 
the dynamics of delay-coupled neurons and investigate the 
period of oscillations which can be induced for sufficiently 
large delay and coupling strength. The additional appli- 
cation of time-delayed self-feedback leads to complex sce- 
narios of synchronized in-phase or antiphase oscillations, 
bursting patterns, or amplitude death. The stabilization 
of steady states in optical systems is studied in Sec. 2 in 
the framework of a modified Lang-Kobayashi model of a 
semiconductor laser where the time-delayed feedback is 
realized by a Fabry-Perot (FP) resonator. 

1. NEUROSYSTEMS 

In order to grasp the complicated interaction between 
billions of neurons in large neural networks, those are often 
lumped into groups of neural populations each of which 
can be represented as an effective excitable element that 
is mutually coupled to the other elements (Rosenblum and 
Pikovsky, 2004; Popovych et al., 2005). In this sense the 



simplest model which may reveal features of interacting 
neurons consists of two coupled neural oscillators. Each of 
these will be represented by a simplified FitzHugh-Nagumo 
system (FitzHugh, 1961; Nagumo et al., 1962), which is 
often used as a paradigmatic generic model for neurons, 
or more generally, excitable systems (Lindner et al., 2004). 
Here we use two identical FitzHugh-Nagumo systems with 
parameters corresponding to the excitable regime. 

We consider the simultaneous action of delayed coupling 
and delayed self-feedback. Here we choose to apply the 
self-feedback term symmetrically to both activator equa- 
tions, but other feedback schemes are also possible. The 
equations of the system are given by 



eui = ui 



■vi C[u2{t - r) - ui{t)] 



-^K[ui{t - tk) - ui{t)] 



vi = ui -\- a 



eu2 = U2 



V2 ^ C[ui{t - r) - U2{t)] 



-^K[u2{t - Tk) - U2{t)] 



V2=U2+ a. 



(1) 

(2) 



(3) 
(4) 



where subsystems Eqs. (l)-(2) and (3)-(4) represent two 
different neurons (or neuron populations), Ui {i = 1,2) 
describing the activator (e.g., transmembrane voltages) 
and Vi modelling the inhibitor (e.g., electrical conduc- 
tances of the relevant ion currents across the respective 
membranes). Here a is a bifurcation parameter whose 
value defines whether the system is excitable (a > 1) or 
demonstrates periodic firing, i.e., autonomous oscillations 
(a < 1), e is positive parameter that is usually chosen to be 
much smaller than unity, corresponding to fast activator 
variables ui, U2, and slow inhibitor variables vi, V2. 

Let us illustrate the dynamics of a single neuron model by 
considering an uncoupled subsystem Eqs. (1) and (2) (C = 
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Fig. 1. Panel (a): Phase portrait with nuhchne of the sin- 
gle FitzHugh-Nagumo system without time-delayed 
feedback. Panel (b) Time series of the activator u 
and inhibitor v as red and green curves, respectively. 
Parameters: e = 0.005, a = 1.05, and noise D = 0.02. 

0, K = 0) under the influence of noise. This is realized by 
small random fluctations modeled as Gaussian white noise 
D({t) with noise intensity D applied to Eq. (2). We fix 
D = 0.02, and also set a = 1.05, e = 0.005. In Fig. 1(a) 
the dashed blue and solid red curves show the nullclines 
of Eq. (1) and (2), respectively, which intersect at a fixed 
point. The phase point that is initially placed at the fixed 
point stays in its close vicinity if the applied random 
perturbation remains small. However, if the perturbation 
is larger than some threshold value, the phase point makes 
a large excursion in the phase space before returning to the 
vicinity of the fixed point again. In Fig. 1(a) the black solid 
line illustrates such a phase trajectory and in Fig. 1(b) 
realizations of ui and vi time series from Eqs. (1) and (2) 
are shown. The motion of the phase point consists of two 
stages: an activation time during which the system waits 
for a sufficiently large perturbation before it can make an 
excursion, and the excursion itself. The excursion time is 
almost completely defined by the deterministic properties 
of the system and is hardly influenced by noise. 

The synaptic coupling between two neurons in Eqs. (1)- 
(4) is modelled as a diffusive coupling considered for 
simplicity to be symmetric (Liley and Wright, 1994; Pinto 
et al., 2000; De-Miguel et al., 2001). More general delayed 
couplings are considered in (Buric and Todorovic, 2003). 
The coupling strength C summarizes how information is 
distributed between neurons. The mutual delay r in the 
coupling is motivated by the propagation delay of action 
potentials between the two neurons ui and U2. 

Besides the delayed coupling we also consider delayed 
self-feedback in the form suggested by Pyragas (Pyragas, 
1992), where the difference s{t) — s{t — tk) of di system 
variable 5, e.g., activator or inhibitor, at time t and at 
a delayed time t — tk, multiplied by some control am- 
plitude is coupled back into the same system. Such 
feedback loops might arise naturally in neural systems. 



e.g., due to neurovascular couplings that has a charac- 
teristic latency, or due to finite propagation speed along 
cyclic connections within a neuron sub-population, or they 
could be realized by external feedback loops as part of 
a therapeutical measure, as proposed by Popovych et al. 
(2005). This feedback scheme is simple to implement and 
quite robust. One distinct advantage of this method is its 
noninvasiveness, i.e., in the ideal deterministic limit the 
control force vanishes on the target orbit, which may be a 
steady state or a periodic oscillation of period r. In case of 
noisy dynamics the control force, of course, does not vanish 
but still remains small, compared to other common control 
techniques using external periodic signals, for instance, 
in deep-brain stimulation to suppress neural synchrony in 
Parkinson's disease (Tass, 2002). 

By a linear stability analysis, it can be shown that the 
fixed point remains stable for all values of K and tk in 
case of a > 1, as without self- feedback (Scholl et al., 2008). 
Redefining ^ = 1 — a^ — C — i^T, one obtains the factorized 
characteristic equation 

1 - + eA^ = XKe-^^^ ± XCe'^^ (5) 

Substituting the Hopf condition X = ico and separating 
into real and imaginary parts yields for the imaginary part 

K cos(cc;rx) ± C cos(cjr) (6) 
This equation has no solution for a > 1 since |^| = — 



5 
4 

3 i 



10 

8 

2 




amplitude 






' — ' — ' — ' — ' — ' 

(a) 12 3 4 5 


' — 

,(b) 


r 




10 

8 

6 



0.2 0.4 0.6 0.8 1 

C 

Fig. 2. Regime of oscillations in the (C, r) parameter 
plane for initial conditions corresponding to single- 
pulse-excitation in one system. The oscillation period 
T is color coded. The transition between black and 
color marks the bifurcation line. Inset (a) shows the 
oscillation period vs. r in a cut at C = 0.8. Inset 
(b): schematic plot of the saddle- node bifurcation of 
a stable (red solid line) and unstable (blue dashed) 
limit cycle. The maximum oscillation amplitude is 
plotted vs. the delay time r and the stable fixed point 
is plotted as a solid black line. The grey background 
marks the bistable region. Parameters: a = 1.05, 
e = 0.01, K = (Schoh et al., 2008). 

In Fig. 2 the regime of oscillations is shown in the param- 
eter plane of the coupling strength C and coupling delay 
r. The oscillation period is color coded. The boundary of 
this colored region is given by the minimum coupling delay 



Tmin as a function of C. For large coupling strength, Tmin 
is almost independent of C; with decreasing C it sharply 
increases, and at some small minimum C no oscillations 
exist at all. At the boundary, the oscillation sets in with 
finite frequency and amplitude as can be seen in the insets 
of Fig. 2 which show a cut of the parameter plane at 
C = 0.8. The oscillation period increases linearly with r. 
The mechanism that generates the oscillation is a saddle- 
node bifurcation of limit cycles (see inset (b) of Fig. 2), 
creating a pair of a stable and an unstable limit cycle. The 
unstable limit cycle separates the two attractor basins of 
the stable limit cycle and the stable fixed point. 
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Fig. 3. Infiuence of delayed self-feedback upon coupled 
oscillations. The mean interspike interval (ISI) is color 
coded in the control parameter plane of the self- 
feedback gain K and delay tk- White areas mark 
regimes of irregular oscillations where the ISI variance 
becomes large (> 0.01). Time series corresponding to 
points (a)-(f) are shown in Fig. 4. Other parameters: 
a = 1.3, e = 0.01, C = 0.5, r = 3 (Schoh et al., 2008). 

The adopted form of control allows for the synchronization 
of the two cells not only for identical values of r and 
Tx, but generates an intricate pattern of synchronization 
islands or stripes in the (r, r^) control parameter plane 
(Fig. 3) corresponding to single-spike in-phase and an- 
tiphase oscillations with constant interspike intervals, see 
also Fig. 4(a)-(d). Further, for adequately chosen param- 
eter sets of coupling and self-feedback control, we observe 
effects such as bursting patterns Fig. 4(f) and oscillator 
death Fig. 4(e). In addition to these effects, there exists 
a control parameter regime in which the self-feedback has 
no effect on the oscillation periods (shaded yellow). 

Fig. 3 shows the control parameter plane for coupling 
parameters of the uncontrolled system in the oscillatory 
regime {C = 0.5 and r = 3). We observere three principal 
regimes: (i) Control has no effect on the oscillation period 
(yellow), although the form of the stable limit cycle is 
slightly altered (Fig. 4(a)). (ii) Islands of in-phase and 
antiphase synchronization (color coded, see Fig. 4 (b)- 
(d)).(iii) Oscihator death (black) Fig. 4 (e)). 

2. LASERS 

Semiconductor lasers with external optical feedback from 
a mirror can be described by the Lang-Kobayashi (LK) 
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4. Different modes of oscillation corresponding to 
different self-feedback parameters K, r (red solid 
lines: activators Ui{t), green solid lines: inhibitors 
Vi{t)). (a), (b): Antiphase oscillations for (a) K = 
0.05, = 3 (period T = 6) and (b) K = 0.5, = 2 
(T = 2); (c), (d): In-phase oscillations for (c) K = 
0.5, = 3 (period T = 3) and (d) K = 0.5, tk = 1.5 
(T = 1.5); (e): Oscillator death for K = 0.9, = 0.9; 
(f): Bursting pattern for K = 0.5, = 3.2. Other 
parameters: a = 1.3, e = 0.01, C = 0.5, r = 3 (Scholl 
et al., 2008). 

model (Lang and Kobayashi, 1980). In dimensionless form, 
it consists of two differential equations for the slowly 
varying amplitude (envelope) E{t) of the complex electric 
field and the reduced carrier density (inversion) n{t). 

Here we consider a modification of the LK equations 
appropriate for multisection semiconductor lasers with an 
internal passive dispersive refiector (Tronciu et al., 2000). 
This is modeled by a gain function k{n) depending upon 
the internal dispersive feedback from the Bragg grating. 
Such a laser structure allows for more complex dynamic 
behavior including self-sustained relaxation oscillations 
(intensity pulsations) generated by Hopf bifurcations, as 
has been shown in the framework of traveling wave laser 
models (Bauer et al., 2004; Schikora et al., 2006). We 
are interested in the regime above a supercritical Hopf 
bifurcation where the fixed point in the uncontrolled 
system is unstable. Combining the rate equation for the 
carrier density from Ref. (Tronciu et al., 2000) with the 
rate equation for the complex electric field, we obtain the 
following form of modified LK equations: 



dE 

~dt 
dn 

'dt 



-n - (l + n)/c(n) \E[ 



(7) 
(8) 



where a denotes the linewidth enhancement factor, / is the 
reduced excess injection current, T is the time scale ratio 



of the carrier lifetime Tc and the photon hfetime r^, and 
Eiy{t) denotes the feedback term, which wih be described 
in detail later. 

The function /c(n), which models the internal dispersive 
feedback, is chosen as a Lorentzian, as proposed by Tronciu 
et al. (2000): 

4(n - no) + 

where A denotes the height, W is the width, and no is 
the position of the resonance. The parameter ko is chosen 
such that k{0) = 1 at the laser threshold. Througout the 
following we will use the parameters A = 1, W = 0.02, 
and no = -0.034. 
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Fig. 5. Schematic diagram of a semiconductor laser with 
resonant feedback from a Fabry- Perot resonator. K 
denotes an attenuator, R is related to the mirror 
reflectivity of the external resonator, r is the round 
trip time of the resonator, cp and S are phase shift 
and latency time due to the distance between laser 
and resonator (Dahms et al., 2008). 

In Eqs. (7), the feedback term E}){t) has not been specifled 
yet. In the following, we introduce the feedback term such 
that it models a Fabry- Perot resonator. As opposed to the 
original LK model where only a single external mirror is 
considered, we take an external FP resonator with multiple 
reflections into account: 



[E{t -S-mT)-E{t-S-{m^ l)r)] 



Ei,{t) = i^^R'' 

= K [E{t - S) - E{t - S - r)] + REi,{t - r). (10) 

r is the delay time (cavity round trip time), is a memory 
parameter (mirror reflectivity), 6 denotes the latency time 
originating from a single round trip between the laser and 
the resonator (see Fig. 5), and = Ke'^"^ is the complex 
feedback gain with amplitude K and the feedback phase 
if which results from the associated optical phase shift. 
Throughout this work we use resonant feedback from the 
FP resonator. 

The latency time 5^ i.e., the propagation time between 
the laser and the FP, is correlated to the phase by the 
relation (f = Qq^^ where is the frequency of the emitted 
light. However, we consider the two parameters (p and S as 
independent variables because the phase (p can be tuned by 
subwavelength changes of the separation between laser and 
FP resonator, on which scale the slowly varying amplitude 
E, which depends upon ^, does not change. The effect of 
latency in time-delayed feedback was already studied in 
a general context in Refs. (Just et al., 1999; Hovel and 
Socolar, 2003; Hovel and Schoh, 2005; Dahms et al., 2007). 

To investigate the stability of the cw laser emission, we 
will perform a linear stability analysis of the lasing state 
is located at {n = 0^ E = >/7e*^) where the phase of the 



electric fleld can be arbitrarily flxed to ip 
laser mode). 



(solitary 



Using the abbreviations E{t) = VT [Qq -\- x{t) -\- iy{t)], 
no = ^/TjT and r = T-i [1 + / (1 + £ J] with real- 
valued X and the flxed point is located at (n = 0, x = 
O^y = 0). An exponential ansatz exp (At) for all three 
variables and n leads to the characteristic equation 



= (2r + A) 



Ke 



-A6 



1 



-At 



1 - Re-^^ 



+2AKe" 



-A6 



-At 



-A5 



1 



— T— COS UD 

Re-^^ ^ 

-At 
1X7 COS (/9 - 



A^ 



- 4nl (A 



(11) 



aKe 



-A5 



1 



-At 



sm (p 



1-Re-^^ ^ I -Re 

In our simulations, we use the following parameters, which 
were chosen close to the values by Tronciu et al. (2000): 
VLq = 0.06, a = 5, r = -0.01 (corresponding to T = 500, 
/ = 1.8, A = 1,W = 0.02, no = -0.034, and ko = 0.993). 
Thus the intrinsic period of the uncontrolled unstable 
focus is To ~ 7r/r^o ^ 52. 

The stability of the lasing flxed point is given by Eq. (11). 
We solved this equation via Newton's method. Since the 
transcendental equation has an inflnite number of roots, 
we scanned the complex plane as initial conditions of 
the root-flnding algorithm to locate the eigenvalue with 
the largest real part. The parameter space consisting of 
r, R, (p, and 6 is flve-dimensional for flxed F and 
Qq. To visualize the domain of control, we consider two- 
dimensional sections of this parameter space. 
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Fig. 6. Domain of control according to Eq. (11) in the 
{K^ r)-plane for different values of cp. The color code 
denotes the largest real part Re (A) of the eigenvalues 
A; only negative values are plotted. Panels (a), (b), 
(c), and (d) correspond to (p = 0, 7r/16, 7r/8, and 
37r/16, respectively. Other parameters: F = —0.01, 
= 0.06, a = 5, = 0.7, and 5 = (Dahms et al., 
2008). 

In Fig. 6, the domain of control is shown in the (i^, r)- 
plane for different values of the phase: (p = 0, 7r/16, 7r/8, 
and 37r/16 in panels (a), (b), (c), and (d), respectively. The 
color code denotes the largest real part of the eigenvalues 
and is therefore a measure of stability. Note that only 




values of Re (A) < are plotted, thus the shaded regions 
correspond to a stable lasing fixed point, i.e., a stable cw 
output. The control domains form tongues separated by 
(white) regions of no control around r = uTq with n 
integer, just as in the generic model studied in Refs. (Hovel 
and Scholl, 2005; Dahms et al., 2007). It can be seen 
that the domain of control shrinks with increasing phase. 
Here, the domains of control are cut off by boundaries 
from the upper left and right for increasing phase, leading 
to overall smaller regions of stability. The tongues of 
stabilization are also slightly distorted towards smaller 
values of r. Additionally, in this picture, the regions of 
optimum stability, denoted by yellow color, are shifted 
towards larger values of the feedback gain K. 
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Fig. 7. Domain of control in the (i^, r)-plane for different 
values of 6 and fixed (p = 0. The color code denotes 
the largest real part Re (A) of the eigenvalues A; only 
negative values are plotted. Panels (a), (b), (c), and 
(d) correspond to 5 = 0, 2.5, 5, and 7.5, respectively. 
Other parameters as in Fig. 6 (Dahms et al., 2008). 

Next, we will investigate the role of the latency time in 
the {K^ r)-plane. In Fig. 7, the domain of control in the 
(i^, r) -plane is depicted for different values of the latency 
time, i.e., ^ = 0, 2.5, 5, and 7.5, and fixed (p = 0. For 
larger latency times, the domains of control shrink and 
it can also be observed that they are bent down towards 
smaller values of the time delay r. Note that the regions 
of optimum stability, denoted by yellow color, are only 
slightly affected by the change of the latency time. 

All figures shown here were obtained for a fixed value of 
the memory parameter = 0.7 as used by Schikora et al. 

(2006) . To investigate the dependence of the control on 
we display the domains of control in the (i^T, i?)-plane 

for different values of the phase {(f = 0, 7r/8, 7r/4, and 
Stt/S) and fixed time delay r = 26 in Fig. 8. This value 
of r was chosen based on the results from the generic 
model considered by Hovel and Scholl (2005); Dahms et al. 

(2007) , where it was shown that the optimum time delay 
is given by r = To/2 = 7r/Im(Ao), where Aq denotes the 
eigenvalue of the uncontrolled system. In the LK model, 
the imaginary part of the eigenvalues in the uncontrolled 
system is given by Eq. (11) with K = 0, i.e., Im(Ao) = 



which yields for our parameters Topt ^ 26. 



Y^4r^Q — r^. This leads to an optimum time delay: 
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Fig. 8. Domain of control in the (i^, i?)-plane for different 
values of if and fixed optimum time delay r = 26. The 
color code denotes the largest real part Re (A) of the 
eigenvalues A; only negative values are plotted. Panels 
(a), (b), (c), and (d) correspond to = 0, 7r/16, 7r/8, 
and 37r/16, respectively. Other parameters as in Fig. 6 
(Dahms et al., 2008). 

Now, in Fig. 8, it can be seen that the domain of control 
in the (i^, i?)-plane has maximum size for (/p = for 
this choice of the time delay r. [See panel Fig. 8(a).] 
For increasing phase, the domain of control shrinks while 
moving to the upper right. Stability is then only achieved 
in a small region at large values of K and R. [See panels 
Fig. 8(b) to (d).] 
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Fig. 9. Domain of control in the (K, (^)-plane for different 
values of r and fixed ^ = 0. The greyscale (color code) 
denotes the largest real part Re (A) of the eigenvalues 
A; only negative values are plotted. Panels (a) and 
(b) correspond to r = 26 and 5, respectively. Other 
parameters as in Fig. 6 (Dahms et al., 2008). 

To investigate the dependence of the domain of control 
on the choice of the phase ip further, we consider another 
two-dimensional projection of the five-dimensional control- 
parameter space parameterized by feedback gain K and 
the feedback phase if. This section is depicted in Fig. 9 
for two different values of the time delay and fixed (5 = 0. 
In panel (a), the time delay is chosen as 26, which is the 
optimum r according to Eq. (12). Here, it can be seen that 
the optimum phase is located at slightly negative values 
for small values of the feedback gain up to K ^ 0.05. 
Increasing the optimum phase changes its sign and is 
now located at small positive values of cp. For the case 
of r = 5, which is depicted in panel (b), stability is 



overall enhanced drastically. The yellow areas, denoting 
regions of optimum stability, are located at negative 
for small K up to K ^ 0.1. Control is possible even for 
a small value of if below — 7r/2, if the feedback gain is 
tuned exactly to the small range of K ^ 0.01. For larger 
feedback gain with K > 0.1, the optimum value of (f is 
located at positive values. The region of optimum stability 
is located at large values of K around 0.2. The shape of 
the control domain in Fig. 9 is markedly different from 
that in the generic normal form model (see Figs. 6(c) and 
7(c) in Ref. (Dahms et al., 2007)), but appears to be in 
line with full device simulations within a travelling wave 
model (Wiinsche et al., 2008). 

To summarize, we have focused on simple rate equation 
models of neural systems and lasers, which contain a time 
delay. In these systems, analytical results are important to 
complement computer simulations to determine whether 
these systems fall into categories defined by the nature 
of their bifurcations. This way, we have gained a more 
general understanding of the principles of dynamical pat- 
terns and their control. Delay dynamics can, on one hand, 
add to the dynamic complexity and generate bifurcations, 
and, on the other hand, stabilize unstable states and 
suppress complex dynamics. Even in the most generic 
neural systems described by the FitzHugh-Nagumo equa- 
tion, when coupled in a minimum network of two nodes 
with delayed coupling and delayed self-feedback, we found 
complex scenarios of synchronized neuronal dynamics with 
in-phase or antiphase oscillations, bursting patterns, and 
amplitude death. In the case of the multisection semicon- 
ductor laser with an internal passive dispersive reflector 
and external optical feedback described by the modified 
Lang-Kobayashi model, we found that time-delayed self- 
feedback can stabilize the laser. In the regime above a 
supercritical Hopf bifurcation, where the fixed point in 
the uncontrolled system is unstable, unwanted intensity 
pulsations can be suppressed by time-delayed feedback 
control. 
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